Introduction to Digital Soili Mapping

My objective is to: (describe specific aim of project and the target soil classes) You may choose to explore your data in Steps 2 and 3 below before deciding on a final objective

My objective for this project is to use random forest to predictively model soil classes for the Essex County, VT data set. I will explorer the data provided, and create additional data as needed. I will then evaluate statistical data reduction techniques by evaluating model performance in an iterative approach to try and maximize model accuracy and minimize OOB estimates of error. Once a final selection of covariates is made, I will use the selected covariates to employ a cost- constrained cLHS to generate a set of sampling points. Finally, I will make model predictions and analyze the outputs.

• note - throughout this document some lines of code were commented out using #, to reduce the length of the final document. These lines of code are mostly reviewing and inspecting data objects.

# Setting up the r environment
# Load the required libraries
library(sp)
library(maps)
library(rgdal)
library(raster)
library(caret)
library(maptools)
library(doParallel)
library(googleVis)
library(psych)
library(clhs)
library(ggplot2)
library(gridExtra)
library(grid)
library(corrplot)

# Set working directory
setwd("C:/Essex/modeling")

Covariates - Identify the characteristics of the study area that can be represented by spectral, terrain, or other ancillary data. Use the SCORPAN framework to organize important study area characteristics and existing data layers that might be used to represent these characteristics. Depending on your objective, it may not be necessary to consider all SCORPAN factors.

Do you need to create any additional data layers to represent the biophysical characteristics identified in question 3? If so, list the data layer and the biophysical characteristic it represents.

The final step in processing was to set all rasters to the same size (resolution), pixel alignment, and extent. This was accomplished using a batch project function in ArcGIS, with snap raster, mask, extent, projection and cell size set to match the 5m DEM provided. Batch project was utilized because during review of the covariate data, MRVBF was found to have some projection issues, not allowing it to be combined in a raster stack in R.

Based on the SCORPAN factors and data listed above, what covariate set did you choose for modeling? Did you use a particular covariate selection method? List covariates and describe selection method below. You may choose to use your expert knowledge and understanding of soil-landscape relationships in lieu of a quantitative method because those weren’t comprehensively covered in this intro course

Random forest models are fairly robust, thus initial modeling was accomplished by including all of the provided covariates plus the additional covariates created. Then several statistical techniques were used to filter down the covariates to increase model accuracy and decrease the out of bag error rate.

Statistical Covariate Reduction Methods:

Near Zero Variance: This method computes the variance of each covariate and returns those that are very near or equal to zero, thus removing data that does not have a high degree of variability. The selected covariates are then removed from the list.

Correlation Reduction: This method generates a correlation matrix of the covariate set. A threshold is selected representing the degree of correlation allowed between the covariates. Those covariates with correlation higher than the selected threshold are removed from the list.

Recursive Feature Elimination: This method runs multiple random forest models and keeps track of the variable importance during each iteration. The covariates that are repeatedly selected as “most” important are returned.

Model 1: This model serves as a baseline. All of the components are modeled within the Essex_field_points shapefile using all developed covariates.

# load Raster Data as a Raster Stack

# read in raster layers to create raster stack for random forest prediction
r.list=list.files(getwd(), pattern="tif$", full.names = FALSE)
r.list

##  [1] "claymin.tif"     "curv.tif"        "esx_dcs.tif"    
##  [4] "esx_maxcurv.tif" "esx_mincurv.tif" "esx_mpwetsm.tif"
##  [7] "esx_mrvbf.tif"   "esx_plncurv.tif" "esx_procurv.tif"
## [10] "esx_saf.tif"     "esx_slope30.tif" "esx_slope60.tif"
## [13] "esx_totsol.tif"  "esx_tri.tif"     "feox.tif"       
## [16] "gndvi.tif"       "mincomp.tif"     "ndvi.tif"       
## [19] "sagawi.tif"      "spi.tif"         "tancurv.tif"    
## [22] "tpi.tif"

#create raster stack
r.stack <- stack(r.list)

#remove files from working environment keeping only the raster stack and shapefile pedon locations (memory management)

to.remove <- ls()
matches <- c("r.stack", "comp$")
to.remove <- c(to.remove[!grepl(paste0(matches, collapse = "|"), to.remove)], "to.remove")
rm(list=to.remove)
#gc()


# set up dataframe for modeling
# first extract raster values from the raster stack by the points

# extract values from r.stack.all and set up a new dataframe
all.df <- extract(r.stack, comp, df=T, na.rm=TRUE)   

# inspect dataframe
#head(all.df)
#names(all.df)
#summary(all.df)
#describe(all.df) # nice method for looking at stats including skewness from library(psych)

# add component information to the dataframe
all.df$class <- comp@data$FEATURE

# inspect dataframe
#head(all.df)
#names(all.df)

# prepare data for  model training

# drop column ID, and Move colum Class to front
# get number of colums from dataframe
#length(all.df)

# output from above
#> length(all.df)
#[1] 26

# dropping column ID. use the length from above for the first number, subtract 1 for the second
comp.rf <- cbind(all.df[ ,26], data.frame(all.df[ ,2:25]))

# inspect comp
#names(comp.rf)

# rename first column
names(comp.rf)[c(1)] <- c('class')

# inspect comp, class should be the first column, followed by the 25 covariates
#names(comp.rf)
#summary(comp.rf)

# inspect naming convention for components
#unique(comp.rf$class)

# need to convert the component names so that there are no spaces
comp.rf$class <- make.names(comp.rf$class)

#unique(comp.rf$class)

# check the data type of the class field
#data.class(comp.rf$class)

# class is a charactor vector needs to be a factor 
# converting class to a factor
comp.rf$class <- factor(comp.rf$class)

#data.class(comp.rf$class)


# Modeling the soil components by class using the Caret Package


#The following code allows you to set the seed to run RF from the caret pckage in parallel and get reproducable results

#The simulation will fit models with subset sizes of 50, 49...1
subsets <- c(1:50)

# set seeds to get reporducable results when running the process in parallel
set.seed(12)
seeds <- vector(mode = "list", length=76)
for(i in 1:75) seeds[[i]] <- sample.int(1000, length(subsets) + 1)
seeds[[76]] <- sample.int(1000, 1)

# in the Train control you can set up the number of folds and repeats, as well as the data splitting, in this instance I have p set to .7, 70 percent of the data will be used for training and 30 percent used for testing.

# set up the train control 
# this is the only instance of fitControl, since it is used over and over keep in working environment
fitControl <- trainControl(method = "repeatedcv", 
                           number = 15,
                           repeats = 5,
                           p = 0.7, #30% used for test set, 70% used for training set
                           selectionFunction = 'best', 
                           classProbs = T,
                           savePredictions = T, 
                           returnResamp = 'final',
                           search = "random",
                           seeds = seeds)


# Random Forest - Parallel process, highlight and run from c1 to stopCluster(c1)
c1 <- makeCluster(detectCores()-1, type='PSOCK')# utilizes all but 1 core
registerDoParallel(c1)
set.seed(48)
rfFit = train(class ~ ., data = comp.rf,
              "rf", 
              trControl = fitControl, 
              ntree = 500, #number of trees default is 500, which seems to work best anyway. 
              tuneLength=10, 
              metric = 'Kappa', 
              na.action=na.pass,
              keep.forest=TRUE, # added this line and the next for partial dependance plots
              importance=TRUE)
stopCluster(c1)
#gc()


# inspect rfFit
rfFit

## Random Forest 
## 
## 185 samples
##  24 predictor
##   8 classes: 'Brayton', 'Cabot', 'Colonel', 'Dixfield', 'Peacham', 'Skerry', 'Skerry.swpd', 'Wilmington' 
## 
## No pre-processing
## Resampling: Cross-Validated (15 fold, repeated 5 times) 
## Summary of sample sizes: 174, 172, 172, 171, 173, 171, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    1    0.5490343  0.3793141
##    3    0.5601114  0.4015874
##    6    0.5476402  0.3851055
##   10    0.5606323  0.4055248
##   12    0.5473489  0.3875904
##   17    0.5584211  0.4023872
##   19    0.5462313  0.3864348
##   24    0.5525797  0.3944286
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 10.

# accuracy 0.5606323

# Get print out of final model
rfFit$finalModel

## 
## Call:
##  randomForest(x = x, y = y, ntree = 500, mtry = param$mtry, importance = TRUE,      keep.forest = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 10
## 
##         OOB estimate of  error rate: 44.86%
## Confusion matrix:
##             Brayton Cabot Colonel Dixfield Peacham Skerry Skerry.swpd
## Brayton           0     4       4        0       0      0           0
## Cabot             2    14       6        0       0      0           0
## Colonel           0     4      36       15       0      0           0
## Dixfield          0     0      15       36       0      0           0
## Peacham           0     1       0        0       0      0           0
## Skerry            0     0       0        5       0      5           0
## Skerry.swpd       0     0       0        0       0      1           0
## Wilmington        0     1      11        5       0      0           0
##             Wilmington class.error
## Brayton              0   1.0000000
## Cabot                3   0.4400000
## Colonel              5   0.4000000
## Dixfield             1   0.3076923
## Peacham              0   1.0000000
## Skerry               0   0.5000000
## Skerry.swpd          0   1.0000000
## Wilmington          11   0.6071429

#OOB estimate of  error rate: 44.86%

Exploration of the data revealed water features within the modeling area. Random forest will try to fit the soil components to the water features. To reduce class error and increase model performance, training points with the class of water were added to the training data set. This was accomplished by editing the Essex_field_points shapefile in ArcGIS using Landsat8 imagery to visually locate the water features and add training points in those locations. This method works best for those components, typically miscellaneous areas, which are easily mapped using imagery.

Model 2: This model is the same as the previous one, but with water features included.

# load component data including water component
comp2 <- readShapePoints("comp2.shp")

# inspecting shapefile
#head(comp2)
#summary(comp2)
#names(comp2)

# look up coordinate system
#comp2@proj4string

# prj2epsg.org

# Assign projection
proj4string(comp2) <- CRS("+init=epsg:32145")

# inspect proojection
#comp2@proj4string

# set up data frame for modeling
# first extract raster values from the raster stack by the points

# extract values from r.stack.all and set up a new dataframe
all.df2 <- extract(r.stack, comp2, df=T, na.rm=TRUE)   

# inspect dataframe
#head(all.df2)
#names(all.df2)
#summary(all.df2)
#describe(all.df2) # nice method for looking at stats including skewness

# Add component information to the dataframe
all.df2$class <- comp2@data$FEATURE

# inspect dataframe
#head(all.df2)
#names(all.df2)

# prepare data for model training

# drop column ID, and Move colum class to front
# get number of colums from dataframe
#length(all.df2)

# output from above
#> length(all.df2)
#[1] 26

# dropping column ID. use the length from above for the first number, subtract 1 for the second
comp.rf2 <- cbind(all.df2[ ,26], data.frame(all.df2[ ,2:25]))

# inspect comp
#names(comp.rf2)

# rename first column
names(comp.rf2)[c(1)] <- c('class')

# inspect comp, class should be the first column, followed by the 25 covariates
#names(comp.rf2)
#summary(comp.rf2)

# inspect naming convention for components
#unique(comp.rf2$class)

# need to convert the component names so that there are no spaces.
comp.rf2$class <- make.names(comp.rf2$class)
#unique(comp.rf2$class)

# check the data type of the class field
#data.class(comp.rf2$class)

# class is a charactor vector needs to be a factor 
# converting class to a factor
comp.rf2$class <- factor(comp.rf2$class)

#data.class(comp.rf2$class)


# Modeling the soil components by class using the Caret Package

# Random Forest - Parallel process, highlight and run from c1 to stopCluster(c1)
c1 <- makeCluster(detectCores()-1, type='PSOCK')# utilizes all but 1 core
registerDoParallel(c1)
set.seed(48)
rfFit2 = train(class ~ ., data = comp.rf2,
              "rf", 
              trControl = fitControl, 
              ntree = 500, #number of trees default is 500, which seems to work best anyway. 
              tuneLength=10, 
              metric = 'Kappa', 
              na.action=na.pass,
              keep.forest=TRUE, # added this line and the next for partial dependance plots
              importance=TRUE)
stopCluster(c1)
#gc()


# Inspect rfFit
rfFit2

## Random Forest 
## 
## 194 samples
##  24 predictor
##   9 classes: 'Brayton', 'Cabot', 'Colonel', 'Dixfield', 'Peacham', 'Skerry', 'Skerry.swpd', 'water', 'Wilmington' 
## 
## No pre-processing
## Resampling: Cross-Validated (15 fold, repeated 5 times) 
## Summary of sample sizes: 183, 181, 180, 179, 183, 181, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    1    0.5803058  0.4393162
##    3    0.5797543  0.4429485
##    6    0.5710730  0.4345155
##   10    0.5795912  0.4461765
##   12    0.5827957  0.4508411
##   17    0.5818753  0.4478714
##   19    0.5756472  0.4406079
##   24    0.5712695  0.4351036
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 12.

# accuracy 0.5827957 

# Get print out of final model
rfFit2$finalModel

## 
## Call:
##  randomForest(x = x, y = y, ntree = 500, mtry = param$mtry, importance = TRUE,      keep.forest = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 12
## 
##         OOB estimate of  error rate: 44.33%
## Confusion matrix:
##             Brayton Cabot Colonel Dixfield Peacham Skerry Skerry.swpd
## Brayton           0     4       4        0       0      0           0
## Cabot             2    13       6        0       1      0           0
## Colonel           0     4      33       14       0      0           0
## Dixfield          0     0      14       37       0      1           0
## Peacham           0     1       0        0       0      0           0
## Skerry            0     0       0        5       0      5           0
## Skerry.swpd       0     0       0        0       0      1           0
## water             0     0       0        0       0      0           0
## Wilmington        0     2      10        5       0      0           0
##             water Wilmington class.error
## Brayton         0          0   1.0000000
## Cabot           0          3   0.4800000
## Colonel         1          8   0.4500000
## Dixfield        0          0   0.2884615
## Peacham         0          0   1.0000000
## Skerry          0          0   0.5000000
## Skerry.swpd     0          0   1.0000000
## water           9          0   0.0000000
## Wilmington      0         11   0.6071429

#OOB estimate of  error rate: 44.33%

Model performance increased slightly with the addition of the water components. The next step is to optimize model performance with the implementation of covariate reduction.

Model 3: This model is the same as the previous one, with the addition of covariate reduction by statistical means. The two methods that follow include near zero variance filtering and correlation reduction.

# convert raster stack to a data frame
stack.df <- as.data.frame(r.stack, na.rm=TRUE)

length(stack.df) # the number of covariates we have in the raster stack

## [1] 24

# 1: Near Zero Variance Filtering

# find covaraiates with near zero variance
# 
# removes covariates whose variance is near or equal to zero
# this process is parallelized to enable faster processing
cl <- makeCluster(detectCores()-1) # makes a cluster using all but 1 cpu cores
registerDoParallel(cl)
# the actual zeroVar function, creates a vector containing which covariates should be removed
zeroVar <- nearZeroVar(stack.df, foreach = TRUE, allowParallel = TRUE)
stopCluster(cl)

# removing unused items from memory
#gc()

# check the zeroVar object
head(zeroVar)

## [1] 2

# *note if integer(0), means that there are no covariates to be removed

# remove covariates with zeroVar
stack.df <- if(length(zeroVar) > 0){
  stack.df[, -zeroVar]
} else {
  stack.df
}

length(stack.df) # the new number of covariates, zeroVar removed 1 covariate

## [1] 23

# nearZero removed one covariate
names(stack.df) # curvature was removed

##  [1] "claymin"     "esx_dcs"     "esx_maxcurv" "esx_mincurv" "esx_mpwetsm"
##  [6] "esx_mrvbf"   "esx_plncurv" "esx_procurv" "esx_saf"     "esx_slope30"
## [11] "esx_slope60" "esx_totsol"  "esx_tri"     "feox"        "gndvi"      
## [16] "mincomp.1"   "mincomp.2"   "mincomp.3"   "ndvi"        "sagawi"     
## [21] "spi"         "tancurv"     "tpi"

# 2. Correlation Reduction

# find correlations within covariates
# a correlation matrix is determined and covariates with high degree of correlation are returned

# create correlation matrix
cor.mat <- cor(stack.df)

# visually examin the correlation matrix
#print(cor.mat)
corrplot(cor.mat)

# many of the curverature covariates and slope derivatives are highly correlated

# I don't want to remove too many variables so I will set the threshold for highly correlated covariates high

# find high degree of correlation, cutoff is the threshold to set. If cutoff = 0.8 then covariates > or = 80 correlated are removed
highCorr <- findCorrelation(cor.mat, cutoff = 0.95)

# remove covariates with high degree of correlation
stack.df <- if(length(highCorr) > 0){
  stack.df[, -highCorr]
} else {
  stack.df
}

length(stack.df) # number of covariates left, corr reduction removed 3 covariates

## [1] 20

names(stack.df)

##  [1] "esx_dcs"     "esx_maxcurv" "esx_mincurv" "esx_mpwetsm" "esx_mrvbf"  
##  [6] "esx_plncurv" "esx_procurv" "esx_saf"     "esx_slope60" "esx_totsol" 
## [11] "esx_tri"     "feox"        "gndvi"       "mincomp.1"   "mincomp.2"  
## [16] "ndvi"        "sagawi"      "spi"         "tancurv"     "tpi"

# Model with remaining covairates

# Try running the model again with the the 20 covariates remaining
# set up data frame for modeling

# create a subset of the original raster stack
r.stack.2 <- subset(r.stack, names(stack.df))

names(r.stack.2)

##  [1] "esx_dcs"     "esx_maxcurv" "esx_mincurv" "esx_mpwetsm" "esx_mrvbf"  
##  [6] "esx_plncurv" "esx_procurv" "esx_saf"     "esx_slope60" "esx_totsol" 
## [11] "esx_tri"     "feox"        "gndvi"       "mincomp.1"   "mincomp.2"  
## [16] "ndvi"        "sagawi"      "spi"         "tancurv"     "tpi"

# extract values from r.stack.2 and set up a new dataframe
all.df3 <- extract(r.stack.2, comp2, df=T, na.rm=TRUE)

# inspect dataframe
#head(all.df3)
#names(all.df3)
#describe(all.df3)
#summary(all.df3)

# add component information to the dataframe
all.df3$class <- comp2@data$FEATURE

# inspect dataframe
#head(all.df3)
#names(all.df3)

# prepare data for model training

# drop column ID, and Move colum class to front
# get number of colums from dataframe
#length(all.df3)

# output from above
#> length(all.df.2)
#[1] 22

# dropping column ID. use the length from above for the first number, subtract 1 for the second
comp.rf3 <- cbind(all.df3[ ,22], data.frame(all.df3[ ,2:21]))

# inspect comp
#names(comp.rf3)

# add class as the first column
names(comp.rf3)[c(1)] <- c('class')

# inspect comp, Class should be the first column, followed by the 43 covariates
#names(comp.rf3)
#describe(comp.rf3)   
#summary(comp.rf3)

# need to convert the component names so that there are no spaces.

comp.rf3$class <- make.names(comp.rf3$class)

#data.class(comp.rf3$class)

# class is a charactor vector needs to be a factor 
# converting class to a factor
comp.rf3$class <- factor(comp.rf3$class)
#data.class(comp.rf3$class)


# Modeling the soil components by class using the Caret Package

# Random Forest - Parallel process, highlight and run from c1 to stopCluster(c1)
c1 <- makeCluster(detectCores()-1, type='PSOCK')
registerDoParallel(c1)
set.seed(48)
rfFit3 = train(class ~ ., data = comp.rf3,
              "rf", 
              trControl = fitControl, 
              ntree = 500, #number of trees default is 500, which seems to work best anyway. 
              tuneLength=10, 
              metric = 'Kappa', 
              na.action=na.pass,
              keep.forest=TRUE, # added this line and the next for partial dependance plots
              importance=TRUE)
stopCluster(c1)
#gc()


# 5.1.2 Inspect rfFit
rfFit3

## Random Forest 
## 
## 194 samples
##  20 predictor
##   9 classes: 'Brayton', 'Cabot', 'Colonel', 'Dixfield', 'Peacham', 'Skerry', 'Skerry.swpd', 'water', 'Wilmington' 
## 
## No pre-processing
## Resampling: Cross-Validated (15 fold, repeated 5 times) 
## Summary of sample sizes: 183, 181, 180, 179, 183, 181, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    1    0.5879715  0.4486188
##    3    0.5890804  0.4529288
##    5    0.5855977  0.4504398
##    8    0.5918594  0.4606135
##   10    0.5776705  0.4423673
##   14    0.5827847  0.4499521
##   16    0.5832954  0.4501917
##   20    0.5862783  0.4542458
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 8.

# accuracy 0.5918594 

# Get print out of final model
rfFit3$finalModel

## 
## Call:
##  randomForest(x = x, y = y, ntree = 500, mtry = param$mtry, importance = TRUE,      keep.forest = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 8
## 
##         OOB estimate of  error rate: 42.27%
## Confusion matrix:
##             Brayton Cabot Colonel Dixfield Peacham Skerry Skerry.swpd
## Brayton           0     4       4        0       0      0           0
## Cabot             2    15       5        0       0      0           0
## Colonel           0     2      36       13       0      0           0
## Dixfield          0     0      15       36       0      0           0
## Peacham           0     1       0        0       0      0           0
## Skerry            0     0       0        5       0      5           0
## Skerry.swpd       0     0       0        0       0      1           0
## water             0     0       0        0       0      0           0
## Wilmington        0     2      11        4       0      0           0
##             water Wilmington class.error
## Brayton         0          0   1.0000000
## Cabot           0          3   0.4000000
## Colonel         0          9   0.4000000
## Dixfield        0          1   0.3076923
## Peacham         0          0   1.0000000
## Skerry          0          0   0.5000000
## Skerry.swpd     0          0   1.0000000
## water           9          0   0.0000000
## Wilmington      0         11   0.6071429

#OOB estimate of  error rate: 42.27%

In model 3, the accuracy increased and OOB error rate decreased. There are several other methods for filtering down covariates. Next, Recursive Feature Elimination will be used to further reduce the covariate set.

Model 4: This model is the same as above, but with Recursive Feature Elimination (RFE) implemented.

#  If data include levels that have only 1 cases,  it does not work... 
#  Gives  Error in { : task 3 failed - "Can't have empty classes in y."
#  So, for RFE we will only use those classes that have more than one case per class

#  Select levels that have more than 1 cases to proceed with the rfe 
subset(table(comp.rf3$class), table(comp.rf3$class) > 1)

## 
##    Brayton      Cabot    Colonel   Dixfield     Skerry      water 
##          8         25         60         52         10          9 
## Wilmington 
##         28

names(subset(table(comp.rf3$class), table(comp.rf3$class) > 1))

## [1] "Brayton"    "Cabot"      "Colonel"    "Dixfield"   "Skerry"    
## [6] "water"      "Wilmington"

# Use this to get names of factors that have more than one case
cat(paste(shQuote(names(subset(table(comp.rf3$class), table(comp.rf3$class) > 1)), type="cmd"), collapse=", "))

## "Brayton", "Cabot", "Colonel", "Dixfield", "Skerry", "water", "Wilmington"

# manually paste the result here and add the 'c(' and the ')' to bookend the list
# manually paste the new list here for subsetting the dataframe
comp.sub <- subset(comp.rf3, class %in% c("Brayton", "Cabot", "Colonel", "Dixfield", "Skerry", "water", "Wilmington"))

nrow(comp.sub)

## [1] 192

# CRITICAL STEP! You must assign the new subset as a factor for this to work
comp.sub$class <- factor(comp.sub$class)

#data.class(comp.sub$class)

Run Recursive Feature Elimination. This selects the covariates that will build the best RF model

# This process is setup to run as a parallel process. Set you number of cpu cores
# in the make cluster function. This example uses detect cores to use all available. To set
# the number of cores just type the number in place of detect cores(). Next, change the subsets


#The simulation will fit models with subset sizes of 43, 42, 41.....1. 
subsets <- c(1:50)

# set seeds to get reporducable results when running the process in parallel
set.seed(12)
seeds <- vector(mode = "list", length=76)
for(i in 1:75) seeds[[i]] <- sample.int(1000, length(subsets) + 1)
seeds[[76]] <- sample.int(1000, 1)

# set up the rfe control
ctrl.RFE <- rfeControl(functions = rfFuncs,
                       method = "repeatedcv",
                       number = 15,
                       repeats = 5,
                       seeds = seeds, 
                       verbose = FALSE)

## highlight and run everything from c1 to stopCluster(c1) to run RFE
#gc()
c1 <- makeCluster(detectCores()-1, type='PSOCK')
registerDoParallel(c1)
set.seed(9)
rf.RFE <- rfe(x = comp.sub[,-1],
              y = comp.sub$class,
              sizes = subsets,
              rfeControl = ctrl.RFE,
              allowParallel = TRUE
)
stopCluster(c1)              
#gc()

# Look at the results
rf.RFE

## 
## Recursive feature selection
## 
## Outer resampling method: Cross-Validated (15 fold, repeated 5 times) 
## 
## Resampling performance over subset size:
## 
##  Variables Accuracy  Kappa AccuracySD KappaSD Selected
##          1   0.4469 0.2853     0.1074  0.1436         
##          2   0.5185 0.3699     0.1272  0.1690         
##          3   0.5277 0.3816     0.1266  0.1681         
##          4   0.5475 0.4086     0.1391  0.1818         
##          5   0.5588 0.4221     0.1214  0.1587         
##          6   0.5726 0.4412     0.1154  0.1511         
##          7   0.5822 0.4525     0.1282  0.1690         
##          8   0.5817 0.4511     0.1235  0.1600         
##          9   0.5955 0.4684     0.1276  0.1672         
##         10   0.5875 0.4568     0.1367  0.1814         
##         11   0.5893 0.4584     0.1395  0.1846         
##         12   0.5973 0.4684     0.1420  0.1886         
##         13   0.5940 0.4647     0.1335  0.1757         
##         14   0.6000 0.4707     0.1287  0.1717         
##         15   0.5934 0.4618     0.1244  0.1675         
##         16   0.5942 0.4638     0.1284  0.1696         
##         17   0.6020 0.4736     0.1223  0.1624         
##         18   0.6058 0.4784     0.1197  0.1567         
##         19   0.6092 0.4823     0.1246  0.1661        *
##         20   0.5943 0.4627     0.1174  0.1548         
## 
## The top 5 variables (out of 19):
##    esx_mrvbf, esx_tri, esx_saf, sagawi, esx_slope60

plot(rf.RFE, metric="Accuracy", main='RFE Accuracy')

# See list of predictors
predictors(rf.RFE)

##  [1] "esx_mrvbf"   "esx_tri"     "esx_saf"     "sagawi"      "esx_slope60"
##  [6] "esx_mpwetsm" "spi"         "esx_dcs"     "esx_maxcurv" "tancurv"    
## [11] "gndvi"       "ndvi"        "mincomp.1"   "esx_totsol"  "feox"       
## [16] "esx_plncurv" "tpi"         "mincomp.2"   "esx_procurv"

# get list of covariates from rfe
predictors(rf.RFE)

##  [1] "esx_mrvbf"   "esx_tri"     "esx_saf"     "sagawi"      "esx_slope60"
##  [6] "esx_mpwetsm" "spi"         "esx_dcs"     "esx_maxcurv" "tancurv"    
## [11] "gndvi"       "ndvi"        "mincomp.1"   "esx_totsol"  "feox"       
## [16] "esx_plncurv" "tpi"         "mincomp.2"   "esx_procurv"

# Create a Random Forest model utilizing the covariates selected using RFE


# create new raster stack to for model development from the predictors above

r.stack.rfe <- subset(r.stack.2, predictors(rf.RFE))
names(r.stack.rfe)

##  [1] "esx_mrvbf"   "esx_tri"     "esx_saf"     "sagawi"      "esx_slope60"
##  [6] "esx_mpwetsm" "spi"         "esx_dcs"     "esx_maxcurv" "tancurv"    
## [11] "gndvi"       "ndvi"        "mincomp.1"   "esx_totsol"  "feox"       
## [16] "esx_plncurv" "tpi"         "mincomp.2"   "esx_procurv"

# extract data data from raster stack
# set up data frame for modeling
# first extract raster values from the raster stack by the points

# extract values from r.stack.all and set up a new dataframe
all.df4 <- extract(r.stack.rfe, comp2, df=T, na.rm=TRUE)   

# inspect dataframe
#head(all.df4)
#names(all.df4)
#summary(all.df4)
#describe(all.df4) # nice method for looking at stats including skewness

# add component information to the dataframe
all.df4$class <- comp2@data$FEATURE

# inspect dataframe
#head(all.df4)
#names(all.df4)

# prepare data for feature selection and model training

# drop column ID, and Move colum class to front
# get number of colums from dataframe
#length(all.df4)

# output from above
#> length(all.df4)
#[1] 21

# dropping column ID. use the length from above for the first number, subtract 1 for the second
comp.rf4 <- cbind(all.df4[ ,21], data.frame(all.df4[ ,2:20]))

# inspect comp
#names(comp.rf4)

# rename first column
names(comp.rf4)[c(1)] <- c('class')

# inspect comp, class should be the first column, followed by the 25 covariates
#names(comp.rf4)

# inspect naming convention for components
#unique(comp.rf4$class)

# need to convert the component names so that there are no spaces.
comp.rf4$class <- make.names(comp.rf4$class)
#unique(comp.rf4$class)

# check the data type of the class field
#data.class(comp.rf4$class)

# class is a charactor vector needs to be a factor 
# converting Class to a factor
comp.rf4$class <- factor(comp.rf4$class)
#data.class(comp.rf4$class)


# Modeling the soil components by class using the Caret Package

# The following code allows you to set the seed to run RF from the caret pckage in parallel and get reproducable results

#The simulation will fit models with subset sizes of 50, 49.....1. 
subsets <- c(1:50)

# set seeds to get reporducable results when running the process in parallel
set.seed(12)
seeds <- vector(mode = "list", length=76)
for(i in 1:75) seeds[[i]] <- sample.int(1000, length(subsets) + 1)
seeds[[76]] <- sample.int(1000, 1)

# in the Train control you can set up the number of folds and repeats, as well as the data splitting, in this instance I have p set to .7, 70 percent of the data will be used for training and 30 percent used for testing.

# set up the train control
fitControl <- trainControl(method = "repeatedcv", 
                           number = 15,
                           repeats = 5,
                           p = 0.7, #30% used for test set, 70% used for training set
                           selectionFunction = 'best', 
                           classProbs = T,
                           savePredictions = T, 
                           returnResamp = 'final',
                           search = "random",
                           seeds = seeds)


#Random Forest - Parallel process, highlight and run from c1 to stopCluster(c1)
c1 <- makeCluster(detectCores()-1, type='PSOCK')# utilizes all but 1 core
registerDoParallel(c1)
set.seed(48)
rfFit4 = train(class ~ ., data = comp.rf4,
               "rf", 
               trControl = fitControl, 
               ntree = 500, #number of trees default is 500, which seems to work best anyway. 
               tuneLength=10, 
               metric = 'Kappa', 
               na.action=na.pass,
               keep.forest=TRUE, # added this line and the next for partial dependance plots
               importance=TRUE)
stopCluster(c1)
#gc()


# inspect rfFit
rfFit4

## Random Forest 
## 
## 194 samples
##  19 predictor
##   9 classes: 'Brayton', 'Cabot', 'Colonel', 'Dixfield', 'Peacham', 'Skerry', 'Skerry.swpd', 'water', 'Wilmington' 
## 
## No pre-processing
## Resampling: Cross-Validated (15 fold, repeated 5 times) 
## Summary of sample sizes: 183, 181, 180, 179, 183, 181, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    1    0.6044566  0.4716902
##    3    0.5937680  0.4603635
##    5    0.5895584  0.4566332
##    8    0.5915738  0.4606959
##    9    0.5871754  0.4538010
##   14    0.6006291  0.4724036
##   15    0.5920193  0.4622670
##   19    0.5891234  0.4573186
## 
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 14.

# accuracy 0.6006291 

# get print out of final model
rfFit4$finalModel

## 
## Call:
##  randomForest(x = x, y = y, ntree = 500, mtry = param$mtry, importance = TRUE,      keep.forest = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 14
## 
##         OOB estimate of  error rate: 43.3%
## Confusion matrix:
##             Brayton Cabot Colonel Dixfield Peacham Skerry Skerry.swpd
## Brayton           0     4       4        0       0      0           0
## Cabot             2    14       6        0       0      0           0
## Colonel           1     4      34       13       0      0           0
## Dixfield          0     0      14       37       0      1           0
## Peacham           0     1       0        0       0      0           0
## Skerry            0     0       0        5       0      5           0
## Skerry.swpd       0     0       0        0       0      1           0
## water             0     0       0        0       0      0           0
## Wilmington        0     2      10        5       0      0           0
##             water Wilmington class.error
## Brayton         0          0   1.0000000
## Cabot           0          3   0.4400000
## Colonel         0          8   0.4333333
## Dixfield        0          0   0.2884615
## Peacham         0          0   1.0000000
## Skerry          0          0   0.5000000
## Skerry.swpd     0          0   1.0000000
## water           9          0   0.0000000
## Wilmington      0         11   0.6071429

#OOB estimate of  error rate: 43.3%

Covariate reduction methods such as nearZero variance, correlation reduction, and recursive feature elimination provide useful means to narrowing a list of covariates to those that contribute the most amount of information without being redundant. A combination of these methods helped to reduce class confusion, and increase overall model accuracy. However, this project could use more sample locations, to improve model accuracy. There may be components not yet described. Some of the components with only one location, may need additional sampling, or correlation to another component.